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We address polarization-induced renormalization of molecular levels in solid-state based single- 
molecule transistors and focus on an organic conjugate molecule where a surprisingly large reduction 
of the addition energy has been observed. We have developed a scheme that combines a self- 
consistent solution of a quantum chemical calculation with a realistic description of the screening 
environment. Our results indeed show a large reduction, and we explain this to be a consequence 
of both (a) a reduction of the electrostatic molecular charging energy and (b) polarization induced 
level shifts of the HOMO and LUMO levels. Finally, we calculate the charge stability diagram and 
explain at a qualitative level general features observed experimentally. 



The recent experimental progress in single-molecule 
electronics has resulted in the realization of the 
three- terminal mole cular single-electron transistor 
(SET) 1 ' 2 ' 3 ' 4 ' 5 ' 6 ' 7 ' 8 ' 9 ' 10 ' shown schematically in Figure [TJl. 
The experimental realizations have been based on a 
variety of techniques including junctions made by elec- 
tromigration, mechanical break junctions and cryogenic 
nanogap fabrication. Many indications of the molecule 
being part of the active transport pathway through the 
junction has been observed. An example is the obser- 
vation of the molecular vibrational excitations, which 
serve as a fingerprint for the molecule 5 . There remain, 
however, several unresolved issues in single-molecule 
transport both in the strong coupling limit, where co- 
herent transport theories seems to strongly overestimate 
the current level, as well as in the weak coupling regime, 
where the observed energy gaps are much smaller than 
expected. Experiments on organic molecules have shown 
that the so-called addition energy, which is the difference 
between the molecular ionization potential (IP) and the 
electron affinity (EA), is heavily reduced com pared to 
its gas phase value in single-molecule SETs^HEl 

Reductions of the excitation gaps is well-known from 
other situations. Theoretical studies of semiconduc- 
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FIG. 1: (a) Schematic illustration of a three-terminal single- 
molecule transistor, (b) Energy level alignment in the molec- 
ular junction showing the position of the molecular ionization 
potential (IP) and electron affinity (EA) levels with respect 
to the work functions of the metallic leads. Also the gate elec- 
trode with voltage V g which couples to EA and IP is schemat- 
ically shown. The polarization shifts P+ and P- of the levels 
due to the junction environment are indicated. 



tor/metal interfaces have shown that the band-gap of 
the semiconductor is narrowed near the interface by the 
screening in the metaff^ 2 ! Experiments using photoe- 
mission techniques and scanning tunneling spectroscopy 
to study the electronic structure of single molecules, 
self-assembled monolayers and organic thin films on 
diele ctric and me tal surfaces have shown similar ef- 
fect J 13 ' 14 ' 15 ' 16 ' 17 ' 18 '. The experimental settings of a single- 
molecule SET is to some degree analogous to the situa- 
tion in electrochemical measurements where the equiv- 
alent of the addition energy, the electrochemical gap, is 
well-known to depend on the dielectric properties of the 
surrounding media 19 . However, the screening environ- 
ments are rather different for the electrochemical setup 
and the single-molecule transistor geometry, with one be- 
ing in ionic solutions or organic solvents and the other in 
solid state low temperature environment. A direct com- 
parison is therefore not possible in general. 

It has been suggested that the reduction of the ad- 
dition energy seen in single-molecule SETs is cause d by 
polarization/image charges in the metallic electrode d 3 * 20 *, 
giving rise to a localization of the charges near the metal- 
lic electrodes. Theoretically only a few other studies have 
addressed the polarization induced renormalization of the 
molecular levels in solid state environments and its im- 
plica tions for the electron transport in molecular junc- 
tions 21 ! 22 ! 23 ! 24 !, and the situation is still very much de- 
bated. Therefore, a more realistic and quantitative theo- 
retical description of the surprisingly large effect is called 
for. The purpose of the present letter is to fill out this 
gap and study the influence of the junction environment 
on the the positions of the molecular levels in a realistic 
single-molecule SET. We have developed a scheme that 
includes the polarizable environment in a quantum chem- 
ical calculation, in which the static polarization response 
of the environment and the molecular charge distribution 
is determined self-consistently. Our calculations on the 
conjugated organic molecule used in experiments 3 5 show 
that a large part of the reduction of the addition energy 
can be accounted for by polarization of the environment. 
By using a simplified expression for the addition energy, 
the reduced addition energy can be understood in terms 
screening of the charging energy of the molecule and a 



2 



closing of the HOMO-LUMO gap. 

For a single-molecule SET operating in the Coulomb 
blockade regime, i.e. with a weak tunnel coupling be- 
tween the molecule and the source/drain electrodes, se- 
quential tunneling is the dominating transport mecha- 
nism. In this regime charge transfer through the molecule 
is possible when either the IP or the EA is positioned 
within the bias window. If on the other hand no levels 
are present in the bias window, current is blocked and the 
molecule remain in a fixed charged state. To reach the 
regime where tranport is possible one can either shift the 
IP and EA levels with the gate-voltage or apply a suffi- 
ciently large source-drain bias. This results in a so-called 
charge stability diagram which maps out the molecular 
charged states as a function of source-drain and gate volt- 
age (see e.g. Figure [I]). The addition energy U can be 
extracted from this diagram by measuring the height of 
the central diamond. Since the ionization potential and 
the electron affinity are given by the difference in total 
energy between the neutral molecule (with TV electrons) 
and the cation (TV — 1) and anion (TV + 1), respectively, 

IP = E N ~ X — E N and EA = E N - E N+1 (1) 

the addition energy can be expressed as, 

U = IP - EA = E N+1 + E N - X - 2E N . (2) 

When the molecule is placed in a nanoj unction, charg- 
ing of the molecule induces polarization charge in the 
junction environment. The formation of the polarization 
charge is associated with stabilizing polarization energies 
P + (added hole) and P_ (added electron) for the cation 
and anion, respectively 14 , which shifts the IP and EA rel- 
ative to their gas phase values as illustrated in Figure ^p. 
The resulting reduction of the addition energy is given by 
the sum P = P + + P_ , i.e. U = IP(g) - EA(g) - P. Nat- 
urally, the polarization energy P depends on the screen- 
ing properties and response times of the environment. In 
single-molecule SETs where the typical current level is on 
the order of I ~ pA — nA, the polarization response of 
the metallic electrodes and gate dielectric (given by the 
plasmon frequency ~ 10 15 s _1 and the phonon frequency 
~ 10 13 s _1 respectively) is orders of magnitudes faster 
than the tunneling rate T = I/e ~ 10 7 — 10 10 s _1 for 
electrons, implying that the polarization energy is given 
by the full static response of the environemnt. 

In the Appendix we have provided a general frame- 
work for evaluating total energies of nanoscale systems in 
the presence of a polarizable environment. The main as- 
sumption of our approach is that the polarizable environ- 
ment responds instantaneously to changes in the charge 
state of the molecule, which according to the above con- 
sideration is a reasonable assumption for single-molecule 
SETs. An electrostatic treatment of the environment 
hence suffices and we derive the following effective Hamil- 
tonian for the nanoj unction, 

H = H s + Hp i + H ext (3) 



where H s is the Hamiltonian of a general nanoscale sys- 
tem (in our case a molecule), 

Hpol = jdr p s (r)<S> ind (r) ~ \ J dr (Ps (r))$md(r) (4) 

describes the interaction between the molecular charge 
distribution p s and the polarization charge through the 
induced potential $j n d, and 

H ext = Jdr p s (r)$ ext (r) (5) 

accounts for external voltages applied to the gate, source 
and drain electrodes. The external potential <& ex t satifies 
Laplace equation with boundary-conditions given by the 
applied voltages on the electrodes. The induced potential 
<f>i n d can be obtained via a solution to Poisson's equation 

-V-[e r (r)V*(r)]=47rp fl (r), (6) 

for the potential & = & s + $i n£ j, where e r is the dielec- 
tric constant of the enviroment and <fr s the potential from 
the molecular charge distribution. The present approach 
thus allows for a continuum description of the environ- 
ment combined with a quantum chemical treatment (e.g. 
DFT or Hartree-Fock) of the molecule. In order to ac- 
count for the molecular charge redistribution due to the 
polarization response of the environment the induced po- 
tential must be included in the usual self-consistent cycle 
of e.g. DFT calculations. 

In the present work a semi-empirical method has been 
combined with a finite element treatment of Poisson's 
equation (see Appendix A for details). We note that 
the addition energy we calculate for the isolated OPV5 
molecule (see below), is underestimated with 1.2 eV as 
compared to the DFT value using the B3LYP exchange- 
correlation functional. In spite of this fact, we still expect 
the polarization energies to be accurate, since the inter- 
action with the polarization charge is treated correctly 
in our approach. That this is indeed the case has been 
confirmed by comparison with other methods^. 

We apply here our method to a single-molecule SETs 
based on the thiol-terminated OPV5 molecule, which is 
a organic conjugated molecule consisting of alternating 
phenylene and vinylene groups (see Figure^). In exper- 
imental realizations of OPV5-based SETs both heavily 
reduced addition energies, acces to several redox states, 
molecular vibrational excitations and Kondo effect have 
been observecPEU. The simulated OPV5-SET is illus- 
trated in Figure |2^i with the molecule lying flat on the 
gate dielectric between the source and drain electrodes, 
which are separated by a 3.2 nm gap. The gate elec- 
trode is separated from the molecule by a 5 nm thick 
layer of gate oxide with dielectric constant e r — 10, cor- 
responding to the high- ft dielectric AI2O3 often used in 
experiment. The relatively high dielectric constant of 
AI2O3 ensures a reasonable capacitive coupling to the 
gate electrode. The gold electrodes are modelled by in- 
finitely high metal blocks. The molecule is placed at a 
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FIG. 2: (a) Illustration of the simulated 0PV5-SET with the molecule lying flat on the AI2O3 gate oxide between the source 
and drain electrodes, (b) Molecular structure of OPV5 and isosurface plots (blue: negative, red: positive) of the HOMO and 
LUMO orbitals for the isolated molecule. (c),(d) Polarization response of the nanojunction illustrated by contour plots of the 
induced potential from the polarization charge at the dielectric and electrode interfaces for the neutral molecule (c) and for the 
anion (d). Due to its partially positive charged thiol groups the overall neutral OPV5 molecule induces a negative potential 
in the nearby electrodes and dielectric that reduces the HOMO-LUMO gap (see text). The additional electron of the anion 
results in a significant polarization of the gate dielectric which reduces the charging energy of the molecule. 



distance of 1 A from the surfaces of the source/drain elec- 
trode and the gate oxide. Since the electrostatic image 
plane of a tomic surfaces is located outside the atomic sur- 
face plan d 25 * 2 ^ 1 , this effectively corresponds to a distance 
between the molecule and the surface atoms on the order 
of van der Waals distance (~ 3 A). 

Table [T] summarizes our findings for the addition en- 
ergy, single-particle HOMO-LUMO gap and polarization 
energy in the following three environments: (i) gas phase, 
(ii) as in Figure [2 1 and (iii) molecule placed in the gap 
between two parallel metal surfaces separated by 3.2 nm. 
The polarization energies due to the presence of the junc- 
tion environments results in significant reductions of the 
addition energies relative to their gas phase values. 

For further analysis we shall use the following simpli- 
fied interpretation of the addition energy: starting with 
two neutral molecules then U is the energy cost of trans- 
ferring an electron from one molecule to other (see equa- 
tion ([5])). Since this process involves the promotion of an 
electron from the HOMO in one of the molecules to the 
LUMO in the other molecule, it is suggestive to write the 
addition energy as the HOMO-LUMO gap of the neutral 
molecule, A#^, plus two times the Coulomb energy, E c , 
required to charge a molecule 



U = A 



HL 



2E C . 



(7) 



This is similiar to the expression for the addition energy 
in the constant-interaction model, which has been used 
successfully for conventional quantum dot SETs^. 

In a naive first guess one would expect the reduction of 
U to be mainly a consequence of screening of the charging 
energy E c . However, Table [I] shows that also the HOMO- 
LUMO gaps are reduced in the polarizable environments. 
The origin of this reduction is illustrated in Figure 
Due to the positively charged thiol groups of the overall 
neutral OPV5, a negative electrostatic potential is in- 
duced in the nearby electrodes and dielectric. This com- 
bined with the localization of the HOMO on the thiol 
groups, see Figure [2)3, shifts the HOMO level to higher 
energy. Similar reasoning for the negatively charged car- 
bon backbone and the LUMO leads to a lowering of the 
LUMO level and hence a closing of the HOMO-LUMO 
gap. The charging energy obtained from equation ^ is 
for OPV5 in gas phase E c = 1.08 eV. The screening re- 
sponse of the nanojunction, which is shown for the OPV5 
anion in Figure |2]i, reduces this value to E c = 75 meV. 
We can therefore understand the reduction of the addi- 
tion energy as a consequence of two parallel effects: i) 
a reduction of the HOMO-LUMO gap and ii) screening 
of the Coulomb repulsion on the molecule which lowers 
the charging energy. Since the majority of the reduction 
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environment 


U A H l P 


gas phase 
SET 
gap 


3.27 1.12 - 
0.68 0.53 2.59 
2.08 0.92 1.19 



i.e. 



TABLE I: Calculated addition energies, U, single-particle 
HOMO-LUMO gaps, Ahl, and polarization energies, P, (all 
in eV) for the thiol-terminated OPV5-molecule in the three 
geometries: gas phase (isolated molecule), SET (geometry as 
in Figure 2a), and gap (molecule placed in the gap between 
two infinite parallel metal surfaces). 



is due to the latter effect which is purely electrostatic in 
nature, the reduction for other molecules of the same size 
of OPV5 will be comparable. For smaller molecules we 
have found that the closer proximity of the polarization 
charge enhances the screening of the charging energy, re- 
sulting in a larger absolute reduction. 

The important role of the gate oxide in the reduction 
of U is clearly demonstrated by the large difference in 
the polarization energy between the SET and gap en- 
vironment in Table [TJ More than half of the polariza- 
tion energy of 2.59 eV in the SET environment is due 
to the gate oxide. As the molecule is lying flat on the 
oxide which has almost metallic-like screening proper- 
ties (the image charge strength of a dielectric surface is 
Qs r = £ r /(£r + 1)) this should come as no surprise. Note, 
however, that the polarization energy in the SET envi- 
ronment is highly dependent on the dielectric constant of 
the gate oxide. With a Si02 oxide layer (s r — 3.9) the 
polarization energy is 2.11 eV. Due to its large distance 
to the molecule the gate electrode has no effect on the 
polarization energy. 

As mentioned in the introductory part, the molecu- 
lar levels can be probed by shifting them with the gate 
voltage or opening the source-drain bias window. Single- 
molecule SETs therefore provide a usefull tool for mea- 
suring the energy differences between the molecular levels 
through the charge stability diagram, albeit in an unnat- 
ural environment. We have calculated the charge stabil- 
ity diagram for the OPV5 device by evaluating total en- 
ergies of the neutral, singly charged and doubly charged 
molecule as a function of gate and source-drain voltage. 
The resulting stability diagram shown in Figure [4] is in 
qualitative agreement with experimental results 3 . It is 
characterized by two small diamonds enclosing a big cen- 
tral diamond. The hight of the central diamond is seen 
to be rsj 0.50 eV instead of 0.68 eV as we found for the 
zero bias value of the addition energy in Table [TJ This 
is due to the fact that the HOMO level moves down- 
wards with the applied source-drain voltage, and hence 
decreasing the threshold for pulling out an electron from 
the HOMO. The non-linear edges on the left side of the 
central diamond is a result of this effect. The heights 
of the two smaller diamonds correspond to the addition 
energies of the anion and cation of the OPV5 molecule, 
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0.12 and a LUMO ~ 0.25 



The small addition energies associated with these states 
stem from their half filled frontier orbitals. There- 
fore, when adding/removing an electron to/from the an- 
ion/cation only the charging energy in equation ^ con- 
tributes. The resulting charging energies are ~ 50 meV 
and ~ 85 meV, respectively, showing that due to the 
different spatial distributions of the HOMO and LUMO, 
the charging energy of the cation and anion are not equal, 
which is implicitely assumed in equation 

One important ingredient in understanding the stabil- 
ity diagram is the gate coupling, a = dE mo i/ 'dV g: i.e. 
how much the energy landscape on the molecule changes 
when a voltage is applied to the gate electrode. For 
usual quantum dot devices this is characterized by a sin- 
gle number, which assumes that all states couple equally 
to the gate. For the OPV5-SET considered here this is 
not the case. As shown in Figure [3^, the gate potential 
varies significantly over the extend of the molecule due 
to screening in the metallic electrodes, which results in 
a higher gate coupling to the LUMO compared to the 
HOMO. In the stability diagram this is reflected in the 
different slopes of the diamond edges, which are given 
by the gate couplings to the different charged states of 
the molecule. This has also been observed in a recent 
experiments The slopes of the diamond edges agree well 
with the calculated gate couplings in Figure [3}}, where 
we read off the values « HOMO 
for an oxide thickness of 5 nm 

In conclusion, by using a method where a continuum 
description of the polarizable junction environment is 
combined with a quantum chemical calculation for the 
molecule, we have studied the effect of polarization in 
an OPV5 single-molecule SET. Our results show a sig- 
nificant modification of the addition energy caused by 
both a closing of the HOMO-LUMO gap and screening 
of the intra- molecular Coulomb interactions. Since the 
majority of the reduction is due to the latter effect, the 
reduction obtained here for OPV5, should be general for 
other molecules of same size. From the calculated charge 
stability diagram we explain at a qualitative level the 
origin of alternating diamond sizes, state dependent gate 
couplings and non-linear diamond edges, which all have 
been observed experimentally. 

Our calculations explain a large part of the reductions 
observed exerimentally, but certainly not all. Other ef- 
fects not accounted for in the present work that can re- 
duce the addtion energy even further could be i) a geom- 
etry where the molecule is closer the the metallic elec- 
trodes as compared to our idealized setup ii) polaron for- 
mation upon charging of the molecule, i.e. relaxation of 
the nuclear configuration, which is associated with relax- 
ation energies on the order of ~ 200 — 300 meV 28 iii) a 
correlation induced localization of the added charge near 
the metallic electrodes beyond what can be captured by 
a mean-field approach. 
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FIG. 3: (a) Contour plot of the gate potential for an oxide thickness of 2.5 nm and with 1 V applied to the gate electrode. The 
metallic electrodes, which are held at V, screen the gate potential significantly which results in a stronger coupling to the 
central part of the molecule, (b) Gate couplings as a function of oxide thickness. The gate couplings have been calculated as 
follows: (Mean) by averaging the potential over the atomic positions of the molecule; (HOMO) and (LUMO) by calculating how 
much the levels move; both with 1 V applied to the gate electrode. The localization of the HOMO on the thiol groups results 
in a much lower gate coupling compared to the LUMO, which is delocalized over the carbon backbone of the molecule. For an 
oxide thickness of 5 nm the gate coupling (Mean) of ~ 0.2 is in good agreement with the value reported in the experiment!^. 
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FIG. 4: Charge stability diagram for the OPV5-SET. The 
color indicates the number of levels positioned in the bias 
window (black: 0, red: 1, yellow: 2) and is hence an indirect 
measure of the current level for a given gate and source-drain 
voltage. The Fermi levels of the gold electrodes have been 
placed in the gap of the molecule as illustrated in Figure [T]d. 
In each of the black diamonds current is blocked leaving the 
molecule in the indicated charge states. The addition energies 
of the charge states of the molecule can be read off from the 
heights of the respective diamonds. 
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APPENDIX A: SUPPORTING INFORMATION 
1. Theoretical Framework 

Here we show how to include the polarizable environ- 
ment of a molecular single-electron transistor (SET) in a 
total energy calculation of the molecule. The main idea is 
the following: when treating the metallic leads and gate 
dielectric in a continuum description they enter only in 
the molecular Hamiltonian as an effective potential. 

In order to address the position of the ionization po- 
tential and electron affinity in the molecular junction, the 
full junction Hamiltonian (i.e. molecule, metallic leads, 
gate dielectric and their mutual interaction) must be con- 
sidered 



(Al) 



The first term is the Hamiltonian of a nanoscale sys- 
tem S (here a molecule). In general, H s is a quantum 
mechanical many-body Hamiltonian that can be treated 
with quantum chemical methods of required accuracy. 
For now H s will not be specified further. The second 
term is the Hamiltonian of the polarizable environment 
(E), which accounts for the energy cost due to build-up 
of the polarization charge in the metallic electrodes and 
the gate dielectric. Since the dynamical polarization re- 
sponse of the environment is orders of magnitude faster 
than the electronic tunneling rates between the molecule 
and leads, an electrostatic treatment of the interaction 
will be sufficient. Therefore we replace the environments 
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by their classical electrostatic energies, which are func- 
tionals of the charge distributions and polarizations in 
the leads and in the dielectric. This approximation ne- 
glects the kinetic energy and the exchange-correlation en- 
ergy associated with the charge build up. 

For the metallic environment where the electrons are 
free to move around the electrostatic energy is 

H m = ^JdrJ dv' Pm(r)V c (r - r')p m (r'), (A2) 

where p m is the charge density of the metal and V c (r) = 
l/|r — r'| is the Coulomb interaction (atomic units are 
used throughout this note if not otherwise stated) . In the 
dielectric environment the polarization charge is bound in 
small dipoles p^, giving rise to a macroscopic polarization 
P = ^ i Pi of the dielectric. The energy of a polarized 
dielectric can be expressed as 

Hd = \ [ d rf d *'pb(r)Vc(r-r')pb(r') 

|P(r)| 2 , (A3) 



X(r,r) 



where p^ = — V • P is the bound charge of the dielectric 
and x_the electric susceptibility. The first term in equa- 
tion (A3) describes the electrostatic dipole-dipole inter- 



action. The second local term, which is equivalent to the 
energy of a spring \kx 2 : , accounts for the energy stored in 

the dipoles p^. Utilizing that V • ^rf^p = S(r — r') this 
term can recast in the form \ fdr j dr f pb(r)V(r; r / )p5(r / ), 
where V is a complicated interaction between the bound 
charges. The total Hamiltonian of the environment can 
now be written 

H E = ±JdvJ dv' p E (v)V(v;v')p E (v'), (A4) 

where p E = p m + pd and V is either V c or V c + V. 
As we shall see in the following, a further specific ation 
of V is not required. The last term in equation (Al) 



accounts for the interaction between the molecule and 
the environment. Since we are focusing on the Coulomb 
blockade regime where the molecule has small tunnel cou- 
plings to the leads, the hybridization term of H SE can be 
neglected. What remains is the Coulomb interaction be- 
tween the spatially separated charges of the molecule and 
of the environment 
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^(rJl^r-r'KM* (A5) 



Notice that the molecular charge has contributions from 
both the ionic cores and the valence electrons of the 
molecule: p s (r) = /w(r) + Pe(r) = Ei ^<K r r i) - 
^t(r)^(r). 

With the Hamiltonian in place we now proceed to elim- 
inate the environment degrees of freedom. Since the part 
of the Hamiltonian involving p E is classical and has no 



dynamics, the solution for p E can be found by minimiz- 
ing with respect p E1 i.e. by setting SH/Sp E = 0. The 
resulting equation can be solved for p E giving (in matrix 
notation) 
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which represents the charge density induced by the sys- 
tem charge p s . Inserting this expression for the induced 
charge of the environment back into the full Hamiltonian 
in equation (Al), the terms involving p E can be recast in 
the form 



H E + H SE = -\ Ps V SE (V^)" 1 V EsPs , (A7) 



in which the degrees of freedoms of the environment 
has been substituted by an effective interaction between 
molecular charges. By inserting the expression in equa- 
tion JA6|) for the induced charge density into equa- 



tion (|A7|), we find that H E + H SE = \p s V SE pi n d, which 



is one half of H SE in equation (A5). This result, which is 



a generalization of the classical image charge problem (a 
point charge placed at a distance z from a perfect con- 
ducting surface where H E + H SE = — l/4z), states that 
the energy cost associated with the build up of the po- 
larization charge is always one half of the energy gained 
by the system S through its interaction with the polar- 
ization charges. 

Due to the presence of the electronic field ope rators 
in p s the effective interaction in equation (A7) must 
be approximated. In our scheme we use a Hartree ap- 
proximation. This is justified because i) correlation ef- 
fects due to the effective interaction are small given the 
absence of a short range interaction and ii) exchange 
is not relevant since self-interactions are possible via 
the image charges. Introducing the induced potential 
$md(r) = fdr'p ind (r')/\r - r'|, the Hartree version of 
equation ( |A7[ ) can be written as 

H e+ H se = J dr Ps(r)®ind(r) ~lj dr (PsWl^W, 

(A8) 

where the last term subtracts the double counted contri- 
butions to the total energy in the first term. 

The induced potential can be obtained by solving Pois- 
son's equation 



V-[ £r (r)V$ tot (r)]=47rp s (r) 



(A9) 



with Dirichlet boundary-conditions on the electrode sur- 
faces Si, i.e. <&tot = if r E . Here, the total potential 
is the sum of the potential from the system charges plus 
the potential from the induced charges: & tot = ^ s -\-^ in d. 

The present approach thus allows us to combine a con- 
tinuum description of the junction environment with a 
quantum chemical description of the molecule. 
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2. Semi-empirical Method 

The molecular part of the total junction Hamiltonian 
in equation ( Al ) is split up into two parts: i) one that ac- 
counts for the isolated molecule and ii) one that accounts 
for polarization/redistribution of charge due to interac- 
tions with the junction environment and/or charging of 
the molecule, 



H s = Ho 



pol • 



(A10) 



In the following our semi-empirical treatment is descriped 
in detail. It is similiar to the one presented in RefP^, 
however, here applied to total energy calculations and 
generalized to include dielectrics in the determination of 
the induced potential 

The part describing the isolated molecule is an effective 
tight-binding Hamiltonian 



^3 



(All) 



where the sums run over atomic valence orbit als {fa} and 
Hs a collective index referring to atom, orbital and spin 
index: i — » ilia. We use Extended Hiickel parameters by 
HoffmanrP^ for the onsite and hopping energies 



ei = -Vi 

tij = — kSij (si + £j ) . 



(A12) 
(A13) 



Here Vi is associated with the ionization energy of the 
valence orbital fa, k is a fitting parameter usually set 
to 1.75 and Sij = (fa\fa) is the overlap between the 
non-orthogonal atomic orbit als. Notice that electron- 
electron interactions are implicitly included in Ho due 
to its parametrized form. 

In the part of the Hamiltonian that accounts for polar- 
ization and charging of the molecule, electron-electrons 
interaction are treated at the Hartree level. Since the 
Hartree potential of the isolated molecule is indirectly ac- 
counted for in Ho, only changes in the Hartree potential 
due to variations in the electron density from its value, 
no, in the isolated molecule are considered 



<SVk(r) = J 



dt 



, 8n(r>) 



(A14) 



where Sn = n — no. Since the Hartree potential depends 
on the electron density, this posses a self-consistent prob- 
lem that must be iterated to convergence. 



To simplify the numerics the integral in equation (A 14) 



is approximated by a sum over atomic point charges given 
by the Mulliken populations n M = Tr \pS] , where p is the 
density matrix 



(A15) 



To avoid problems with the diverging point charge po- 
tential when evaluated at the atomic positions, r M , the 



onsite contribution to the sum is replaced by a species de- 
pendent Hubbard U taking into account the energy cost 
of adding an electron to the atom. These param eters are 
taken from the quantum chemical CNDO metho d 31 * 32 !. In 
order to keep consistency between the onsite and the off- 
site interactions, the Magata-NishimotcP 3 ^ interpolation 
formula is used for the latter 



1 



(A16) 



With these approximation the final form of the Hartree 
potential in equation (A14) becomes 



<5V H (i>) = Sn^ + SnvU^. (A17) 

In our atomic basis the polarization/charging part of the 
Hamiltonian is written 



(A18) 



where the last term substracts double counting in the first 
term and Hq. The matrix representation of the Hartree 
potential has been approximated as follows 



Vi = {4>i |tfWi> « 

V ij = {4> i \8Vn\<P j )^\s ij {V i + V j ). 



(A19) 



The part of the Hamiltonian involving the induced po- 
tential is written similarly 

H E + H es =Y1 V ^3 C l C 3 ~\^2( Z »- U v) ^ind(r^) 

(A20) 

with Vjj — (fa\Vi n d\<f>j) evaluated as above in equa- 
tion flA19|). 



3. Poisson's Equation and &i n d 

In the following section it is described how the in- 
duce d po tential is determined by solving Poisson's equa- 
tion (A9) for the total potential <f>tot- 



a. Finite element approach 

One of the major advantages of the finite element 
method (FEM) is its partitioning of the solution domain 
into a finite number of elements (typically triangles in 
2D and tetrahedra in 3D). The possibility to refine the 
element size around sharp corners and in the vicinity of 
spatially rapidly varying source terms, allows FEM to 
handle a large variety of problems and solution domains 
of practically any geometry, hence making it very suit- 
able for modelling of nanoscale devices. 
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(b) 




FIG. 5: Finite element mesh and simplified junction 
geometry, (a) Surface mesh for the 0PV5-SET. (b) Simpli- 
fied junction for which Poisson's equation can be solved ana- 
lytically. In this geometry the electrodes are modelled by infi- 
nite metallic surfaces and the gate oxide as a semi-infinite di- 
electric layer sandwiched between the electrode surfaces. The 
gate electrode is not included (see text). 



In the present case we seek to solve Poisson's equa- 
tion in the geometry of a single- molecule SET. To this 
end we use the finite element software from the FEniCS 
project 34 . 

The molecular charge density, which has both elec- 
tronic and ionic contributions, is represented by a sum of 
atomic centered gaussian charge distributions: 



; (r) = (27r)- 3 / 2 V^% e -l r - r Ml 2 /2o 



(A21) 



Here a — <2q is the width of the gaussians and = 
%n — (Z^ being the atomic valence) is the net atomic 
charge. 



Having obtained the total potential, &tot, the induced 
potential can be extracted by substracting the potential, 
<£ s , from the molecular source charges, which for a gaus- 
sian charge distribution centered at the origin is: 



<D s (r)— erff 



\V2cr 



(A22) 



where erf is the error function. 

The calculated induced potential has been converged 
(to 0.05 eV) with respect to the element size and the 
spatial dimensions of the device. In order to get an accu- 
rate description of the potential at the atomic positions, 
the molecule is enclosed in a box with a fine mesh of el- 
ement size ~ ao. The element size at the boundaries of 
the device are 20 times larger, which is illustrated by the 
boundary mesh in figure [5| a). The following device size 
was used to converge the potential: an electrode hight of 
50 a.u., a device width of 100 a.u. and a device length of 
150 a.u.. 



b. Analytical solution 

In the simplified device geometry shown in Fig. |5|b), 
Poisson's equation can be solved analytically. It can be 
shown that the Greens function satisfying: 



- V- [e r (r)VG(r,r')] = S(r - r' 
has the following image charge solution: 



(A23) 



a=±l r=±l 



e r + r 
e r + 1 



1 



y/(x - crx') 2 + (y - ry') 2 + (z - z') 2 
1 

^ V yfJ^aL^ (x - crx')) 2 + (y - ry') 2 + (* -7f 
1 

^(2nL + (x- ax')) 2 + (y - ry') 2 + (z - z')< 



E 



(A24) 



where L is the electrode spacing (z is the direction per- 
pendicular to the paper plane in Fig.[5jb)) and the sums 
run over all repeated images of the source charge in the 
metallic surfaces and the dielectric. By leaving out the 
contribution from the source charge itself, only the in- 
duced potential remains: 

$ ind (r)=G(r,r') (A25) 



In the calculations using the analytical solution, the 
molecular charge distribution has been approximated by 
point charges located at the atomic positions. The 
induced potential follows directly from the Greens func- 
tions of the individual point charges: 



Here G denotes the Greens function which does not in- 
clude the (7 = 1 and r = 1 term outside the n-sum in 
equation (A24). 



^W = E«fG(r,r,) (A26) 
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Miolecule 


FEM 


Analytic 


0PV2 


3.56 


3.63 


0PV3 


3.11 


3.17 


0PV4 


2.81 


2.87 


0PV5 


2.59 


2.63 



TABLE II: Calculated polarization energies P (in eV) for 
the two junction geometries illustrated in figure [5] The same 
electrode spacing has been used in the two geometries. Due to 
the infinite metallic surfaces, the polarization energy is slighty 
larger for the simplified junction geometry. 



c. FEM vs analytic 

Table [IT] summarizes the calculated polarization ener- 
gies for the realistic junction modelled with FEM and the 



simplified junction shown in Fig.[5|b). The difference be- 
tween the polarization energies in the two junctions is on 
the order of ~ 50 meV for the different OP V- molecules. 
The slighty larger polarization energies in the simplified 
junction stem from the infinite metallic electrodes, which 
screen the Coulomb interactions on the molecule better 
than the semi-infinite metal blocks of the realistic junc- 
tion. Notice that, due to its large distance to the molecule 
(~ 5 nm) and the metallic-like screening properties of the 
AI2O3 gate oxide, the gate electrode does not contribute 
to the polarization energy in the realistic junction. This 
explains the relative small differences between the polar- 
ization energies in the two junction geometries. 

Hence, our analytical solution provides a realistic 
description of the potential in generic SET geometry 
(Fig. 2(a) in the main part of the paper) that can be 
used instead of computationally heavy Poisson solvers. 
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